    # rm(list=ls()) ; gc()
    
    
    #Article:    The Politicization of Bureaucrats: Evidence from Brazil			
    #Authors:     Anderson Frey <anderson.frey@rochester.edu>	& Rogerio Santarrosa <rogeriobs2@insper.edu.br>		    
    
    # This script produces all Tables and Figures in the Article
    # Tables are produced in a way that they can be directly copied/pasted in a LaTex document
    
    # Produced and replicable with RStudio v2024.04.1 and R v4.2.1

    # It uses the following files:  data_main.csv (Tables 2,3,4,5,7,8 and Figures 3,4)
                                    #data_cadunico.csv (Tables 1,4,5,6,8 and Figure 4)                                    
                                    #data_partisanship.csv (Figure 2)
                                    #data_bureau.csv (Tables 1,8)

    
    # Please upload all libraries and functions below before attempting to produce any Tables or Figures
  
    # The seed for all bootstraps in set to the University of Rochester ZIP CODE: 14627

    # Set the working directory to the folder that contains the replication files, use the function below
 
  
      setwd("")
      
    
    # Install Packages ==============================================================================================================
    
      install.packages("stringr") 
      install.packages("lubridate") 
      install.packages("data.table") 
      install.packages("dplyr") 
      install.packages("stringi")
      install.packages("lfe")
      install.packages("rdrobust")
      install.packages("fastDummies")
      install.packages("ggplot2")
      install.packages("boot")
      install.packages("gridExtra")
    
    # ==============================================================================================================================
    
    
    # Libraries ====================================================================================================================
      
      library(stringr) 
      library(lubridate) 
      library(data.table) 
      library(dplyr) 
      library(stringi)
      library(lfe)
      library(rdrobust)
      library(fastDummies)
      library(ggplot2)
      library(boot)
      library(gridExtra)
    
    # ==============================================================================================================================
    
     
    # Functions ===================================================================================================================

      .LatexTable <- function(data,rownames,spaces,dist){
        cols <- ncol(data)*2+1
        res <- data.frame(matrix(0,nrow(data),cols))
        res[,] <- "&"
        fill <- seq(from=2,to=cols,by=2)
        res[,fill] <- data
        res <- cbind(rownames,res)
        lastone <- rep("\tabularnewline",nrow(res))
        lastone[spaces] <- paste0("\tabularnewline[",dist,"cm]")
        lastone[nrow(res)] <- paste0("\tabularnewline[",dist,"cm]")
        res[,ncol(res)] <- lastone
        names(res) <- NULL
        return(res)
      }
      

      .NewPlot <- function(data,r,v,xvar,minD,maxD,bins1,bins2,shift,bb,x.title,y.title,p.title,degree,multi,oneside,bb3){
        
        
        Di <- matrix(0,(bins1),3)
        Di[,1] <- seq(from=minD, to=maxD, length.out=(bins1))
        label <- paste("t",xvar,sep=":")
        
        for (i in 1:nrow(Di)){
          Di[i,2] <- r["t",1] + r[label,1]*Di[i,1] 
          Di[i,3] <-(  v["t","t"] + (Di[i,1]**2)*v[label,label] + 2*Di[i,1]*v["t",label] )^.5
        }
        
        MT <- data.frame(Di) ; names(MT) <- c("D","coeff","sd")
        MT$h <- MT$coeff+1.96*MT$sd
        MT$l <- MT$coeff-1.96*MT$sd  
        
        jazz <- data.frame(matrix(0,(bins2),2))
        xx <- seq(from=min(MT$D), to=(max(MT$D)), length.out=(bins2+1))
        
        data <- data.frame(data)
        data$x <- data[,match(xvar,names(data))]
        
        i <- 2
        
        for (i in 1:nrow(jazz)){
          top <- xx[i+1]
          bot <- xx[i]
          if(is.na(top))  top <- xx[i+1] + 1
          jj <- subset(data, x >= bot & x < top )
          if(oneside==1) jj <- subset(data, x >= bot & x < top & t == 1)
          if(oneside==1) {datax <- subset(data,  t == 1)}else{datax <- data}
          jazz[i,1] <- (xx[i]+xx[i+1])/2
          
          jazz[i,2] <- nrow(jj)/nrow(datax)
        }
        names(jazz) <- c("D","dens")
        
        bb2 <- bb+shift
        
        ggplot(data=MT, aes(x=D, y=coeff, ymin=0))   +    theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          geom_line(data=MT, aes(x=D, y=coeff+shift),colour="gray10",  linewidth=1, alpha=.75) +
          theme(plot.title = element_text(family = "A", size=11))+
          geom_hline(yintercept=(0+shift), lty=3)+
          xlab(x.title) + ylab(y.title) +
          theme(axis.text = element_text(colour="black", size=11, family="A")) +
          theme(axis.title = element_text(colour="black", size=11, family="A")) +
          #theme(plot.title = element_text(colour="black", size=13, family="A")) +
          geom_line(data=MT, aes(x=D, y=h+shift),colour="gray10",  linetype="dashed",linewidth=.75, alpha=.75) +
          geom_line(data=MT, aes(x=D, y=l+shift),colour="gray10",  linetype="dashed",linewidth=0.75, alpha=.75) +
          geom_bar(data=jazz, aes(x=D, y=dens*multi, group = 1), stat="identity", position="dodge",
                   size =0.1,color="gray35",fill="gray85", alpha=0.3) +
          scale_y_continuous(breaks=bb2, labels=bb )+
          scale_x_continuous(breaks=bb3)+
          coord_cartesian(ylim=c(bb2[1],bb2[length(bb2)]))
      }
      
      
      .DiffRDD1 <- function( cities, base, indices, b1a, b1b, v1, v2, con, byvar){
        
        library(dplyr)
        library(data.table)
        library(lfe)
        
        b <- data.table( cities[indices,] )
        bb <- data.table( inner_join(b,base,by="ibge") )
        
        bbx <- filter(bb, abs(run) < b1a )
        if(byvar==1) bbx <- filter(bbx, i==0 )
        bbx$var <- bbx %>% select(all_of(v1))
        if(con==0) reg <- felm(var ~ t*run  |0|0|0 ,  data=bbx) 
        if(con==1) reg <- felm(var ~ t*run +age+gender+public_servant+newcomer+education+past_memb1+past_memb2+past_memb3+pt_coalition+pt_mayor+pt_federal+size  |0|0|0 ,  data=bbx) 
        if(con==2) reg <- felm(var ~ t*run +age+gender+public_servant+newcomer+education+past_memb1+past_memb2+past_memb3+size  |0|0|0 ,  data=bbx) 
        a <- summary(reg)$coefficients["t",1]
        
        bbx <- filter(bb, abs(run) < b1b )
        if(byvar==1) bbx <- filter(bbx, i==1 )
        bbx$var <- bbx %>% select(all_of(v2))
        if(con==0) reg <- felm(var ~t*run |0|0|0 ,  data=bbx) 
        if(con==1) reg <- felm(var ~t*run +age+gender+public_servant+newcomer+education+past_memb1+past_memb2+past_memb3+pt_coalition+pt_mayor+pt_federal+size |0|0|0 ,  data=bbx) 
        if(con==2) reg <- felm(var ~t*run +age+gender+public_servant+newcomer+education+past_memb1+past_memb2+past_memb3+size |0|0|0 ,  data=bbx) 
        b <- summary(reg)$coefficients["t",1]
        
        return(a-b)
        
      }
    
    # ==============================================================================================================================
   
    
    # Basic Setup | Run before starting Tables and Figures =========================================================================
    
      windowsFonts( A=windowsFont("Liberation Sans")) # Font for ggplot2 
      
      party_list1 <- c(15,22,40,12,65,10,20,19,36,13) # Parties in PT's federal base in 2010
      party_list2 <- c(15,22,40,12,65,10,20,19,36)    # Parties in PT's federal base in 2010 (ex-PT)
      
      controls <- c("age", "gender", "public_servant", "education", "newcomer","pt_mayor", 
                    "pt_federal", "pt_coalition", "size", # define covariates
                    "past_memb1", "past_memb2", "past_memb3")
      

      btype <- "cerrd"  # type of algorithm for bandwidth calculation
      qq <- 2           # polynomial for bandwidth calculation
      vtype <- "hc3"    # type of standard error for the RD estimation
      
    # =============================================================================================================================
    
   
    # Necessary Files | Run before starting Tables and Figures=====================================================================
   
      cadunico <- fread( "data_cadunico.csv")
      MainData <- fread( "data_main.csv")
      
     

    # =============================================================================================================================
      
    
    ## T A B L E S  
      
       
    # Table 1 | Summary Statistics: CadUnico Interviewers and Other Bureaucrats ===================================================
        
        bureau <- fread("data_bureau.csv")
      
        summa1 <- cadunico %>% select(ibge,new_party,old_party)
        summa2 <- bureau %>% select(ibge,new_party,old_party)

        # Manually input the total number of partisans in Brazil, aggregated at country level
        quadro <- data.frame(c(12153.479,416.975,1786.411,128352.147))
        quadro[,2] <- quadro[,1]*100/sum(quadro[,1])
        
        # Run the averages
        for (i in 1:2){
          
          if(i==1) xc <- summa1
          if(i==2) xc <- summa2              
          
          
          xc$par <- ifelse(!is.na(xc$old_party),0,3)
          xc$par[!is.na(xc$new_party) & !is.na(xc$old_party)] <- 1
          xc$par[!is.na(xc$new_party) & is.na(xc$old_party)] <- 2
         
          
          
          xc$count <- 1
          agg <- group_by(xc, par) %>% summarise( count=sum(count)/1000)
          agg$pc <- agg$count*100 / sum(agg$count)
          
          if(i==1){aggx <- agg}else{aggx <- cbind(aggx,agg)}
          
        }
        
        # prepare output
        output <- data.frame(matrix(0,4,6))
        output[1:4,1:2] <- aggx[,2:3]
        output[1:4,3:4] <-   aggx[,5:6]
        output[1:4,5:6] <-   quadro[,1:2]
        output$X5 <- output$X5 - output$X3 - output$X1
        output$X6 <- output$X5*100/sum(output$X5)
        
        
        output[1:4,] <- format(round(output[1:4,],1),nsmall=1)
        cols <- ncol(output)*2+1
        resx <- data.frame(matrix(0,nrow(output),cols))
        resx[,] <- "&"
        
        fill <- seq(from=2,to=cols,by=2)
        resx[,fill] <- output
        
        # Print in LaTex format
        resx[,ncol(resx)] <- c("\tabularnewline","\tabularnewline","\tabularnewline","\tabularnewline[0.2cm]")
        lista_colnames <-  c("Old Partisans","Switchers","New Partisans","Non-Partisans")
        resx <- cbind(lista_colnames,resx)
        print(resx, row.names=FALSE, colnames=FALSE)
        
        rm(bureau,agg,aggx,quadro,output,resx,summa1,summa2,xc)
        
    # =============================================================================================================================
     
  
    # Table 2 | Main Results ======================================================================================================
        
       
        w <- MainData
        w$var1 <- w$outcome*100/w$total
        w$var2 <- w$outcome.2*100/w$total.2
        w$var3 <- w$outcome.3*100/w$total.3
        
        # Create Main Table
        res <- data.frame(matrix(0,5,3))
        for (i in 1:3){
          
          if(i == 1) {w$var <- w$var1}
          if(i == 2) {w$var <- w$var2}
          if(i == 3) {w$var <- w$var3}
          
          bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = qq,bwselect = btype, masspoints="off")$bws ; bands
          b1 <- bands[1]
          b2 <- bands[3]
          
          if(i == 1) b1a <- b1
          if(i == 2) b1b <- b1
          if(i == 3) b1c <- b1
      
          reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype)
          reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90 , vce=vtype)
          
          if(i == 1) c1a <- reg$coef[1]
          if(i == 2) c1b <- reg$coef[1]
          if(i == 3) c1c <- reg$coef[1]
          
          
          int <- round(reg$tau_cl[1],3)
          co <- round(reg$coef[1],3)
          se <- round(reg$se[3],3)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          n <- round(reg$N_h[1],0)
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          
          res[1,i] <- paste(format(co,nsmall=3),stars,sep="")
          res[2,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
          res[3,i] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          
          res[4,i] <- format(int,nsmall=3)
          
          res[5,i] <- format(round(b1,2),nsmall=2)
          res[6,i] <- n
          print(i)
          
        }
        
        res1 <- res
       
        # Run the bootstrap for SEs of the Differences
        cities <- w
        cities <- unique(cities, by="ibge")
        cities$x <- 1
        cities <- cities %>% select(ibge,x)
        
        set.seed(14627)
        results <- boot(data=cities, base=w, statistic=.DiffRDD1 , R=500,con=0,byvar=0,
                        b1a=b1a, b1b=b1b,  parallel = "snow", ncpus = 4, v1="var1", v2="var2")
        
        ci = boot.ci(results, type=c("perc","basic"), conf=.95, vce=vtype)  
        ci2 = boot.ci(results, type=c("perc","basic"), conf=.90, vce=vtype) 
        res2 <- data.frame(matrix(0,6,2))
        
        stars <- ""
        if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
        if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
        
        ci1 <- round(ci$perc[4],3)
        ci2 <- round(ci$perc[5],3)
        
        res2[1,1] <- paste(format(round(c1a-c1b,3),nsmall=3),stars,sep="")
        res2[3,1] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
        res2[2,1] <-  paste( "(" ,format(round(sd(results$t),3),nsmall=3), ")",sep="")
        res2[4:6,1] <- "-"
        
        set.seed(14627)
        results <- boot(data=cities, base=w, statistic=.DiffRDD1 , R=500,con=0,byvar=0,
                        b1a=b1a, b1b=b1c,  parallel = "snow", ncpus = 2, v1="var1", v2="var3")
        
        ci = boot.ci(results, type=c("perc","basic"), conf=.95, vce=vtype)  
        ci2 = boot.ci(results, type=c("perc","basic"), conf=.90, vce=vtype) 
        
        stars <- ""
        if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
        if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
        
        ci1 <- round(ci$perc[4],3)
        ci2 <- round(ci$perc[5],3)
        
        res2[1,2] <- paste(format(round(c1a-c1c,3),nsmall=3),stars,sep="")
        res2[3,2] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
        res2[2,2] <-  paste( "(" ,format(round(sd(results$t),3),nsmall=3), ")",sep="")
        res2[4:6,2] <- "-"
        
        # Print Table for LaTex format
        resx <- cbind(res1,res2)
        rw <- c("Incumbency Effect"," ","C.I. (95%)","Baseline","Bandwidth","Municipalities")
        a <- .LatexTable(resx,rw,spaces=c(3),dist=0.2)
        print(a, row.names=FALSE)
        
        rm(a,ci,bands,w,cities,results,resx,res,res1,res2,reg,reg2)
        
    
    # =============================================================================================================================

       
    # Table 3 | Comparison | Matching =============================================================================================
      
      # Select relevant subsample and redefine variables
      
      wx <- filter(MainData, total.m>0)
      wx$v1 <- wx$outcome.m.i*100/wx$total.m
      wx$v2 <- wx$outcome.m.b*100/wx$total.m
      wx$var <- wx$v2 
      
      #Set bands
      bands <- rdbwselect(wx$var, wx$run, c = 0, p = 1, q = qq,bwselect = btype , masspoints="off")$bws ; bands
      b1 <- bands[1]
      b2 <- bands[3]
 
      # Create Main Table
      res <- data.frame(matrix(0,5,3))
      j <- 1
      for (j in 0:1){
        
        #
        if(j==0) wx$var <- wx$v1
        if(j==1) wx$var <- wx$v2
         
        if(j==0) b1a <- b1
        if(j==1) b1b <- b1
        
        covars <- wx %>% select(size,age,gender,newcomer, public_servant,education,past_memb1,past_memb2,past_memb3,pt_coalition,pt_mayor,pt_federal)
        covars <- cbind(covars)
        
        reg <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars,vce=vtype) 
        reg2 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars, level=90, vce=vtype) 
        reg3 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, vce=vtype) 
        
        if(j==0) c1 <- reg$coef[1]
        if(j==1) c2 <- reg$coef[1]
        co <- round(reg$coef[1],3)
        se <- round(reg$se[3],3)
        ci1 <- round(reg$ci[3,1],3)
        ci2 <- round(reg$ci[3,2],3)
        n <- round(reg$N_h[1],0)
        stars <- ""
        if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
        if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
        bas <- round(reg3$tau_cl[1],3)
        if(j==0) bas1 <- bas
        if(j==1) bas2 <- bas
        
        res[1,j+1] <- paste(format(co,nsmall=3),stars,sep="")
        res[3,j+1] <- format(bas,nsmall=3)
        res[2,j+1] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
        dt <- filter(wx, abs(run) < b1 )
        res[4,j+1] <- length(unique(dt$ibge))
        res[5,j+1] <- format(round(b1,2),nsmall=2) 
      }
   
      # Run the bootstrap for SEs of the Differences
      cities <- wx
      cities <- unique(cities, by="ibge")
      cities$x <- 1
      cities <- cities %>% select(ibge,x)
      
      set.seed(14627)
      results <- boot(data=cities, base=wx, statistic=.DiffRDD1 , R=500,con=1,v1="v1",v2="v2",
                      b1a=b1a, b1b=b1b,  parallel = "snow", ncpus = 4, byvar=0)
      
      ci = boot.ci(results, type=c("perc","basic"), conf=.95, vce=vtype)  ;  ci 
      ci2 = boot.ci(results, type=c("perc","basic"), conf=.90, vce=vtype) ;  ci2
      
      stars <- ""
      if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
      if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
      
      ci1 <- round(ci$perc[4],3)
      ci2 <- round(ci$perc[5],3)
      
      res[1,3] <- paste(format(round(c1-c2,3),nsmall=3),stars,sep="")
      res[2,3] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
      res[3,3] <- paste(format(round(bas1-bas2,3),nsmall=3),sep="")
      res[4:5,3] <- "-"
 
      
      # Print Table in LaTex format
      rw <- c("Incumbency Effect","C.I. (95%)","Baseline","Municipalities","Optimal Bandwidth")
      a <- .LatexTable(res,rw,spaces=c(3),dist="0.2")
      print(a, row.names=FALSE)
      
      rm(a,ci,bands,wx,dt,cities,results,res,reg,reg2,reg3,covars)
      
  
    # =============================================================================================================================
      
     
    # Table 4 | Comparison | Senior vs Junior =====================================================================================
 
        # Select relevant subsamples
        w1 <- w2 <- MainData
        w1$i <- 0
        w2$i <- 1
        w <- rbind(w1,w2)
  
        # Create separate outcomes for senior or junior interviewers
        newzy <- cadunico
        newzy$t <- newzy$coalition_side

        newzy$inter_tot <- 1 
        newzy$inter_rel <- ifelse(!is.na(newzy$new_party) & newzy$date_int < newzy$date_new,1,0)
        newzy$i <- ifelse(newzy$date_int <= "2008-09-30"  ,0,1 )
        aa <- data.table(group_by(newzy,ibge,t,i) %>% summarise(outcome.new=sum(inter_rel), inter_tot=sum(inter_tot) ))
        bb <- data.table(aa %>% group_by(ibge,i) %>% summarise(total.new=sum(inter_tot)))
        
        w <- left_join(w,bb,by=c("ibge","i"))
        w <- left_join(w,aa,by=c("ibge","t","i"))
        w[is.na(w)] <- 0
        w$var <- w$outcome.new*100 / w$total.new
        w <- filter(w, !is.na(var))
        rm(newzy,aa,bb)
        
        # Set bands
        bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = 2,bwselect = btype , masspoints="off")$bws 
        b1 <- bands[1] 
        b2 <- bands[3]
        
        # Create Main Table
        res <- data.frame(matrix(0,5,3))
        j <- 1
        for (j in 0:1){
          
          wx <- filter(w, i == j)  
          if(j==0) b1a <- b1
          if(j==1) b1b <- b1
          covars <- wx %>% select(age,gender,newcomer, size,public_servant,education,past_memb1,past_memb2,past_memb3,pt_coalition,pt_mayor,pt_federal)
          covars <- cbind(covars)
          
          reg <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars,vce=vtype) 
          reg2 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars, level=90, vce=vtype) 
          reg3 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, vce=vtype) 
          
          if(j==0) c1 <- reg$coef[1]
          if(j==1) c2 <- reg$coef[1]
          co <- round(reg$coef[1],3)
          se <- round(reg$se[3],3)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          n <- round(reg$N_h[1],0)
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          bas <- round(reg3$tau_cl[1],3)
          if(j==0) bas1 <- bas
          if(j==1) bas2 <- bas
          
          res[1,j+1] <- paste(format(co,nsmall=3),stars,sep="")
          res[2,j+1] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          res[3,j+1] <- paste(format(bas,nsmall=3),sep="")
          
          dt <- filter(wx, abs(run) < b1 )
          res[4,j+1] <- length(unique(dt$ibge))
          res[5,j+1] <- format(round(b1,2),nsmall=2) 
        }
     
        # Run the bootstrap for SEs of the Differences
        cities <- w
        cities <- unique(cities, by="ibge")
        cities$x <- 1
        cities <- cities %>% select(ibge,x)
        
        set.seed(14627)
        results <- boot(data=cities, base=w, statistic=.DiffRDD1 , R=500,con=1,v1="var",v2="var",
                        b1a=b1a, b1b=b1b,  parallel = "snow", ncpus = 4, byvar=1)
        
        ci = boot.ci(results, type=c("perc","basic"), conf=.95, vce=vtype) 
        ci2 = boot.ci(results, type=c("perc","basic"), conf=.90, vce=vtype) 
        
        stars <- ""
        if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
        if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
        
        ci1 <- round(ci$perc[4],3)
        ci2 <- round(ci$perc[5],3)
        
        res[1,3] <- paste(format(round(c1-c2,3),nsmall=3),stars,sep="")
        res[2,3] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
        res[3,3] <- paste(format(round(bas1-bas2,3),nsmall=3),sep="")
        res[4:5,3] <- "-"
        res1 <- res
      
        # Print Table for LaTex format
        rw <- c("Incumbency Effect","C.I. (95%)","Baseline","Municipalities","Optimal Bandwidth")
        a <- .LatexTable(res,rw,spaces=c(2,5),dist="0.2")
        print(a, row.names=FALSE)
        
        rm(a,w,wx,cities,dt,w1,w2,covars,bands,ci,results,res,res1,reg,reg2,reg3)
        
      
    # =============================================================================================================================
    
           
    # Table 5 | Comparison | Number of Interviews ====================================================================================
        
        # Select relevant subsamples
        wx <- MainData
        w1 <- w2 <- wx
        w1$i <- 0
        w2$i <- 1
        wxx <- rbind(w1,w2)
        
        # Get data on the number of interviews per bureaucrat
        newzy <- cadunico
        newzy$t <- newzy$coalition_side
        newzy$inter_tot <- 1 
        newzy$inter_rel <- ifelse(!is.na(newzy$new_party) & newzy$date_int < newzy$date_new,1,0)
        
        # Create the 2 Main Tables
        for (k in 1:2){
          
          
          wxx <- rbind(w1,w2)
          if(k==1) newzy$i <- ifelse(newzy$int_200912 <= newzy$median1  ,1,0 )
          if(k==2) newzy$i <- ifelse(newzy$int_2012 <= newzy$median2  ,1,0 )
          
          aa <- data.table(group_by(newzy,ibge,t,i) %>% summarise(outcome.new=sum(inter_rel),inter_tot=sum(inter_tot)))
          bb <- data.table(group_by(aa,ibge,i) %>% summarise(total.new=sum(inter_tot)))
          aa$inter_tot <- NULL
          
          wxx <- left_join(wxx,bb,by=c("ibge","i"))
          wxx <- left_join(wxx,aa,by=c("ibge","t","i"))
          wxx[is.na(wxx)] <- 0
          wxx$var <- wxx$outcome.new*100 / wxx$total.new
          wxx <- filter(wxx, !is.na(var))
          
          bands <- rdbwselect(wxx$var, wxx$run, c = 0, p = 1, q = 2,bwselect = btype , masspoints="off")$bws 
          b1 <- bands[1] 
          b2 <- bands[3]
          
          
          res <- data.frame(matrix(0,5,3))
       
          for (j in 0:1){
            
            wx <- filter(wxx, i == j)  
      
            if(j==0) b1a <- b1
            if(j==1) b1b <- b1
            
            covars <- wx %>% select(age,gender,newcomer, size,public_servant,education,past_memb1,past_memb2,past_memb3,pt_coalition,pt_mayor,pt_federal)
            covars <- cbind(covars)
            
            reg <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars,vce=vtype) 
            reg2 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars, level=90, vce=vtype) 
            reg3 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, vce=vtype) 
            
            if(j==0) c1 <- reg$coef[1]
            if(j==1) c2 <- reg$coef[1]
            co <- round(reg$coef[1],3)
            se <- round(reg$se[3],3)
            ci1 <- round(reg$ci[3,1],3)
            ci2 <- round(reg$ci[3,2],3)
            n <- round(reg$N_h[1],0)
            stars <- ""
            if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
            if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
            bas <- round(reg3$tau_cl[1],3)
            if(j==0) bas1 <- bas
            if(j==1) bas2 <- bas
            
            res[1,j+1] <- paste(format(co,nsmall=3),stars,sep="")
            res[2,j+1] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
            res[3,j+1] <- paste(format(bas,nsmall=3),sep="")
            
            dt <- filter(wx, abs(run) < b1 )
            res[4,j+1] <- length(unique(dt$ibge))
            res[5,j+1] <- format(round(b1,2),nsmall=2) 
          }
          res
          
          # Run the bootstrap for SEs of the Differences
          cities <- wxx
          cities <- unique(cities, by="ibge")
          cities$x <- 1
          cities <- cities %>% select(ibge,x)
          
          set.seed(14627)
          results <- boot(data=cities, base=wxx, statistic=.DiffRDD1 , R=500,con=1,v1="var",v2="var",
                          b1a=b1a, b1b=b1b,  parallel = "snow", ncpus = 4, byvar=1)
          
          ci = boot.ci(results,  type=c("perc","basic"), conf=.95)  
          ci2 = boot.ci(results, type=c("perc","basic"), conf=.90) 
          
          stars <- ""
          if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
          if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
          
          ci1 <- round(ci$perc[4],3)
          ci2 <- round(ci$perc[5],3)
          
          res[1,3] <- paste(format(round(c1-c2,3),nsmall=3),stars,sep="")
          res[2,3] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          res[3,3] <- paste(format(round(bas1-bas2,3),nsmall=3),sep="")
          res[4:5,3] <- "-"
          if(k==1) res1 <- res
          if(k==2) res2 <- res
        
        }
        
        # Print Panels for LaTex format
        rw <- c("Incumbency Effect","C.I. (95%)","Baseline","Municipalities","Optimal Bandwidth")
        a <- .LatexTable(res1,rw,spaces=c(5),dist="0.2")
        print(a, row.names=FALSE)
        
        rw <- c("Incumbency Effect","C.I. (95%)","Baseline","Municipalities","Optimal Bandwidth")
        a <- .LatexTable(res2,rw,spaces=c(5),dist="0.2")
        print(a, row.names=FALSE)
        
        rm(newzy,wxx,wx,cities,covars,bands,ci,results,res,res1,res2,reg,reg2,reg3,a,aa,bb,dt,w1,w2)

        
    # =============================================================================================================================
        
          
    # Table 6 | Career Upside =====================================================================================================
      
      # Get relevant interviewers  
      xx <- filter(cadunico, coalition_side>=0 & !is.na(remains_bureaucrat) )
 
      xx$ocup12 <- as.numeric(xx$ocup12)
      xx$t <- xx$coalition_side
      xx$ps12 <- xx$remains_bureaucrat
  
      # Create variables
      xx$salary12[is.na(xx$salary12)] <- 0
      xx$v1 <- ((xx$salary)/100)
      xx$v2 <- ((xx$salary12)/100)
      xx$v3 <- xx$v2 / xx$v1 -1
      
      xx$boss1 <- ifelse(xx$ocup < 192000 ,1,0)
      xx$boss2 <- ifelse(xx$ocup12 < 192000 ,1,0)
      xx$boss3 <- xx$boss2 - xx$boss1
 
      xx$i3 <-  xx$civil12 - xx$civil
      
      xxx <- filter(xx,  salary12>0  )
      xd <- filter(xxx,  civil==1 & is.na(old_party)  )
    
      # Create Table
      res <- data.frame(matrix(0,16,3))
      i <- 1; k <- 1
      for (i in 1:4){
        
        xxx$ps <- xxx$ps12
        xx$ps <- xx$ps12
        if(i==1) xx$var <- xx$ps
        if(i==2) xx$var <- xx$ps
        if(i==3) xxx$var <- xxx$v3
        if(i==4) xxx$var <- xxx$boss3
        if(i==1) xd <- xx
        if(i==2) xd <- filter(xx,  civil==0   )
        if(i>=3) xd <- filter(xxx, ps12==1   )
        
        for (k in 1:3){
         if(k==1) xd2 <- xd
         if(k==2) xd2 <- filter(xd,  is.na(old_party) )
         if(k==3) xd2 <- filter(xd,  !is.na(old_party) )
        
        xd2$reg <- substring(xd2$ibge,1,2)
        reg <- felm(var ~ t |  0  | 0 | ibge , data=xd2) ; summary(reg)

        r <- summary(reg)$coefficients
        v <- vcov(reg)
        
 
          coeff <- round(r["t",1],3)
          bas <- round(r["(Intercept)",1],3)
          sd <- round(r["t",2],3)
          p <- (r["t",4])
          st <- ""
          if(p < 0.1) st <- "+"
          if(p < 0.05) st <- "*"
          res[i*4-3,k] <-   paste0(  format( coeff  , nsmall=3) ,st)
          res[i*4-2,k] <-  paste0(  "(",format( sd  , nsmall=3) ,")")
          res[i*4-1,k] <-  paste0(  format( bas  , nsmall=3) )
          
          res[i*4,k] <-  format( round(nrow(xd2),0), nsmall=0)
          
  
        }
        
      }

      res <- res[-c(12),]
 
  
      # Print Table in LaTex format
      rw <- c("Remains in Bureaucracy"," ","Baseline","Observations")#,
      a <- .LatexTable(res[1:4,],rw,spaces=c(4),dist=0.2)
      print(a, row.names=FALSE)
      
      rw <- c("Remains in Bureaucracy"," ","Baseline","Observations")
      a <- .LatexTable(res[5:8,],rw,spaces=c(4),dist=0.2)
      print(a, row.names=FALSE)
      
      rw <- c("Wages"," ","Baseline","Promotion"," ","Baseline","Observations")
      a <- .LatexTable(res[9:15,],rw,spaces=c(3,6,7),dist=0.2)
      print(a, row.names=FALSE)
      
      rm(a,reg,r,v,res,xd2,xx,xd,xxx)
    
    # =============================================================================================================================
        
 
    # Table 7 | Comparison | Parties ==============================================================================================
    

        w <- MainData
        w$var1 <- w$outcome*100/w$total
        
        # Set bands
        bands <- rdbwselect(w$var1, w$run, c = 0, p = 1, q = qq,bwselect = btype , masspoints="off")$bws 
        b1 <- bands[1] 
        b2 <- bands[3]
    
        # Create relevant variables based on partisanship
        w$var <- w$var1
        w$pwin <- w$party*w$t
        w$test <- str_detect(w$colig,"/13/")
        w$pt_mayorb <- ifelse( !w$party %in% party_list1,1,0)
        w$ptcoligwin <- w$ptcolig*w$t
        w <- data.table(w %>% group_by(ibge) %>% mutate(pwin=max(pwin), ptcoligwin=max(ptcoligwin)))
        w$pt_federal_win <- ifelse(w$pwin %in% c(party_list1)  ,1, 0)
        w$pt_mayor_win <- ifelse(w$pwin %in% c(13),1, 0)
        w <- data.table(w %>% group_by(ibge) %>% mutate(pt_there1=sum(pt_mayor), pt_there2=sum(pt_federal), pt_there3=sum(pt_mayorb)))

        # Create Tables
        for (k in 1:2){
     
          if(k==2) {wxx <- filter(w, pt_there2 == 1)  ; wxx$i <- wxx$pt_federal_win }
          if(k==1) {wxx <- filter(w, pt_there1 == 1 ) ; wxx$i <- wxx$pt_mayor_win  } 
  
          res <- data.frame(matrix(0,5,3))
          j <- 1
          for (j in 0:1){
     
            wx <- filter(wxx, i == 1-j)  
            
            if(j==0) b1a <- b1
            if(j==1) b1b <- b1
            
            covars <- wx %>% select(age,gender,newcomer, size,public_servant,education,past_memb1,past_memb2,past_memb3)
            covars <- cbind(covars)
            
            reg <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars,vce=vtype) 
            reg2 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars, level=90, vce=vtype) 
            reg3 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, vce=vtype) 
            
            if(j==0) c1 <- reg$coef[1]
            if(j==1) c2 <- reg$coef[1]
            co <- round(reg$coef[1],3)
            se <- round(reg$se[3],3)
            ci1 <- round(reg$ci[3,1],3)
            ci2 <- round(reg$ci[3,2],3)
            n <- round(reg$N_h[1],0)
            stars <- ""
            if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
            if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
            bas <- round(reg3$tau_cl[1],3)
            if(j==0) bas1 <- bas
            if(j==1) bas2 <- bas
            
            res[1,j+1] <- paste(format(co,nsmall=3),stars,sep="")
            res[2,j+1] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
            res[3,j+1] <- paste(format(bas,nsmall=3),sep="")
            
            dt <- filter(wx, abs(run) < b1 )
            res[4,j+1] <- length(unique(dt$ibge))
            res[5,j+1] <- format(round(b1,2),nsmall=2) 
          }
          res
          
          # Run the bootstrap for SEs of the Differences
          cities <- wxx
          cities <- unique(cities, by="ibge")
          cities$x <- 1
          cities <- cities %>% select(ibge,x)
          
          set.seed(14627)
          results <- boot(data=cities, base=wxx, statistic=.DiffRDD1 , R=500,con=2,v1="var",v2="var",
                          b1a=b1a, b1b=b1b,  parallel = "snow", ncpus = 4, byvar=1)
          
          ci = boot.ci(results, type=c("perc","basic"), conf=.95, vce=vtype)  #;  ci 
          ci2 = boot.ci(results, type=c("perc","basic"), conf=.90, vce=vtype) #;  ci2
          
          stars <- ""
          if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
          if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
          
          ci1 <- round(ci$perc[4],3)
          ci2 <- round(ci$perc[5],3)
          
          res[1,3] <- paste(format(round(c1-c2,3),nsmall=3),stars,sep="")
          res[2,3] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          res[3,3] <- paste(format(round(bas1-bas2,3),nsmall=3),sep="")
          res[4:5,3] <- "-"
          if(k==1) res1 <- res
          if(k==2) res2 <- res
          
        }
        
        # Print Table in LaTex format
        rw <- c("Incumbency Effect","C.I. (95%)","Baseline","Municipalities","Optimal Bandwidth")
        a <- .LatexTable(res1,rw,spaces=c(5),dist="0.2")
        print(a, row.names=FALSE)
        
        rw <- c("Incumbency Effect","C.I. (95%)","Baseline","Municipalities","Optimal Bandwidth")
        a <- .LatexTable(res2,rw,spaces=c(5),dist="0.2")
        print(a, row.names=FALSE)
        
        rm(w,wx,covars,bands,ci,results,res,res1,res2,reg,reg2,reg3,a,wxx,cities,dt)
        
        
        
        
      
        
    # =============================================================================================================================
        
   
    # Table 8 | Candidacy Effects =================================================================================================
  
      bureau <- fread("data_bureau.csv")
      wxx <- MainData
      
      # Get bureaucrats that become candidates
      newzy <- cadunico
      newzy$t <- newzy$coalition_side2
      newzy <- filter(newzy, t >=0)# &
      newzy$inter_rel <- 1
      aa <- data.table(group_by(newzy,ibge,t) %>% summarise(outcome.new=sum(inter_rel)))
      wxx <- left_join(wxx,aa,by=c("ibge","t"))
      wxx[is.na(wxx)] <- 0
      
      newzy <- bureau
      newzy$t <- newzy$coalition_side2
      newzy <- filter(newzy, t >=0)# &
      newzy$inter_rel <- 1
      aa <- data.table(group_by(newzy,ibge,t) %>% summarise(outcome2.new=sum(inter_rel)))
      wxx <- left_join(wxx,aa,by=c("ibge","t"))
      wxx[is.na(wxx)] <- 0
      
      wxx$outcome3.new <- wxx$council_cand - wxx$outcome2.new - wxx$outcome.new
      
      w <- wxx
      rm(newzy,wxx, aa)

      # Create table
      res <- data.frame(matrix(0,6,3))
      for (i in 1:3){
        
        if(i == 1) {w$var <-w$outcome.new*100/w$total}
        if(i == 2) {w$var <- w$outcome2.new*100/w$total.2}
        if(i == 3) {w$var <- w$outcome3.new*100/w$total.3}
        
        bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = qq,bwselect = btype, masspoints="off")$bws ; bands
        b1 <- bands[1]
        b2 <- bands[3]
        
        if(i == 1) b1a <- b1
        if(i == 2) b1b <- b1
        if(i == 3) b1c <- b1
        
        reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype)
        reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90 , vce=vtype)
        
        if(i == 1) c1a <- reg$coef[1]
        if(i == 2) c1b <- reg$coef[1]
        if(i == 3) c1c <- reg$coef[1]
        
        
        int <- round(reg$tau_cl[1],3)
        co <- round(reg$coef[1],3)
        se <- round(reg$se[3],3)
        ci1 <- round(reg$ci[3,1],3)
        ci2 <- round(reg$ci[3,2],3)
        n <- round(reg$N_h[1],0)
        stars <- ""
        if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
        if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
        
        res[1,i] <- paste(format(co,nsmall=3),stars,sep="")
        res[2,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
        res[3,i] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
        
        res[4,i] <- format(int,nsmall=3)
        
        res[5,i] <- format(round(b1,2),nsmall=2)
        res[6,i] <- n
        print(i)
        
      }
      
      # Print Table in LaTex format
      rw <- c("Incumbency Effect"," ","C.I. (95%)","Baseline","Bandwidth","Municipalities")
      a <- .LatexTable(res,rw,spaces=c(6),dist=0.2)
      print(a, row.names=FALSE)
      
      rm(bureau,a,bands,reg,reg2,res,w)
      
    # =============================================================================================================================

      
    ## F I G U R E S 
    

    # Figure 1 | BF Expansion =====================================================================================================
    
      # Number of total BF beneficiaries for the whole country:
      bfb <- c(6571.839,8700.445,10965.810,11043.076,10557.996,12370.915,12946.313,13361.495,13902.155)
      # Target of total BF beneficiaries for the whole country:
      bft <- c(11101.876,11101.876,11101.876,11101.876,11101.876,12995.195,12995.195,13737.202,13737.202)
  
      # Both entered manually, obtained from:
      #Frey, Anderson. 2019. “Cash transfers, clientelism, and political enfranchisement: Evidence from Brazil.”
      #Journal of Public Economics 176:1 – 17.
  
      # Years
      bfy <- 2004:2012
    
      plot_bf <- data.frame(bfb,bft,bfy)
      plot_bf$bft <- plot_bf$bft-5000
      plot_bf$bfb <- plot_bf$bfb-5000
      lab <- c(5000,7500,10000,12500,15000)
      brk <- c(0,2500,5000,7500,10000)
      brk_x <- 2004:2012
    
      ggplot(data=plot_bf, aes(x=bfy, y=bft))   +
        theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
        geom_bar( aes( y=bft),color="gray95", alpha=0.25, stat = "identity")+
        geom_line( aes( y=bfb), color="gray15", alpha=0.75, linewidth=2)+
        scale_y_continuous(labels=lab,breaks=brk)+
        scale_x_continuous(breaks=brk_x) +
        theme(axis.title = element_text(colour="black", size=11, family="A")) +
        theme(axis.title.x = element_blank()) +
        theme(axis.text = element_text(colour="black", size=11, family="A")) +
        ylab("BF Benefitiaries (thousand)")
      
      rm(plot_bf)
    
    # =============================================================================================================================
  
      
    # Figure 2 | Party membership =================================================================================================

      fil_tot <- fread("data_partisanship.csv")
      
      # Aggregate all partisans before 2001
      fil_tot0 <- filter(fil_tot, year<2001)
      fil_tot0 <- group_by(fil_tot0,party) %>% summarise(stock=sum(members) - sum(cancel))
      
      # Total number of voters in teh country by year, entered manually with data from TSE
      a <- sort(rep(c(109782.873,115183.907,121257.350,125603.401,130246.642,135604.038,140394.103),2))
      b <- 2000:2013
      ab <- data.frame(b,a) ; names(ab) <- c("year","total")
      for(i in 2:12){ if(ab$year[i] %% 2 != 0) ab$total[i] <- ab$total[i+1]/2+ab$total[i-1]/2  }
      
      
      # Create Table
      res <- matrix(0,7*12,3)
      
      prt <- c(11,13,15,25,40,45)
      count <- 1
      for (j in 1:6){
        
        
        for (i in 2001:2012){
          xxx <- filter(fil_tot, party == prt[j] & year == i)
          xxx2 <- filter(fil_tot0, party == prt[j] )
          
          res[count,2] <- prt[j]
          res[count,1] <- i
          if(i == 2001) {res[count,3] <- xxx$members - xxx$cancel + xxx2$stock}else
          {res[count,3] <- xxx$members - xxx$cancel + res[(count-1),3]}
          count <- count+1
        }
        
      }
      
      # Plot
      plot_fili <- data.table(res)
      names(plot_fili) <- c("year","party","stock")
      plot_fili <- inner_join(plot_fili, ab, by="year")
      
      plot_fili <- plot_fili %>% arrange(party,year)
      plot_fili$past <- plot_fili$stock
      plot_fili$past[2:nrow(plot_fili)] <- plot_fili$past[1:(nrow(plot_fili)-1)]
      plot_fili$ch <-  (plot_fili$stock - plot_fili$past )/10 / plot_fili$total
      
      pf1 <- filter(plot_fili, year >= 2002 )
      pf1$party <- as.factor(pf1$party)
      
      ggplot(data=pf1, aes(x=year, y=ch, group=party))   +
        theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
        geom_bar( aes(x=year, y=ch, group=party, fill=party), stat = "identity", position = "dodge" ) +
        scale_x_continuous(breaks=2002:2012)+ 
        scale_fill_grey(name  ="Party",labels=c("PP", "PT","MDB","DEM","PSB","PSDB"),start = 0.8, end = 0)+
        theme(axis.title = element_text(colour="black", size=11, family="A")) +
        theme(axis.title.x = element_blank()) +
        theme(axis.text = element_text(colour="black", size=11, family="A")) +
        ylab("New members (% of Brazilian voters)")+
        theme(legend.position=c(0.7,0.65))
      
      rm(ab,fil_tot0,fil_tot,plot_fili,pf1,res,xxx,xxx2,count)
    

    # =============================================================================================================================
      
     
    # Figure 3 | Comparison | CU Expansion ========================================================================================
      
        w <- MainData
        
        # Define variable
        w$var <- w$outcome*100/w$total
        w$tgt3 <- (w$pot_enrollment-mean(w$pot_enrollment))/sd(w$pot_enrollment) 

        # Set bands
        bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = qq,bwselect = btype , masspoints="off")$bws 
        b1 <- bands[1] 
        b2 <- bands[3]
        
        #Run Regression
        data <- filter(w,  abs(run) < b1   )
        data$wght <- b1 - abs(data$run)
        regols <- felm( var ~  t*run*tgt3 +age+gender+public_servant+newcomer+education
                        +past_memb1+past_memb2+past_memb3+pt_coalition+size  | party+reg | 0 | ibge ,weights=data$wght , data=data) 
       
  
        rr <- summary(regols)$coefficients 
        vv <- vcov(regols , robust=T)
    
        vec <- round(seq(from= -12, to = 12,length.out=7),2) 
        vec2 <- round(seq(from= -4.5, to = 3.5,length.out=17),2)
        ph1 <- .NewPlot(data=data,r=rr,v=vv,xvar="tgt3",minD=-4.5,maxD=3.5,bins1=12,bins2=16,shift=12,bb=vec,degree=1,multi=24,oneside=0,
                        x.title="Potential Enrollment (stdevs. from mean)",y.title="Incumbency Effect",p.title="",bb3=vec2) ;   ph1
        
        rm(ph1,rr,vv,bands,w,data,regols)
     
    # =============================================================================================================================
 
  
    # Figure 4 | Comparison | demographic effects==================================================================================
        
     
        w1 <- w2 <- MainData
        w1$i <- 0
        w2$i <- 1
        v <- rbind(w1,w2)
        v$var <- v$outcome*100/v$total
        
        
        # Get individual interviewers
        c <- cadunico 
        c$t <- c$coalition_side
        c <- data.table( c %>% group_by(ibge) %>% mutate(med_salary=median(salary)))
        
        # Set bands
        bands <- rdbwselect(v$var, v$run, c = 0, p = 1, q = qq,bwselect = btype, masspoints="off")$bws 
        b1 <- bands[1]
        b2 <- bands[3]
        
        # Run tables
        res <- data.frame(matrix(0,6,5))
        res2 <- res
        
        for (i in 1:6){
        
          c$inter_tot <- 1 
          c$inter_rel <- ifelse(!is.na(c$new_party) & as.numeric(c$date_int) < (as.numeric(c$date_new))  ,1,0)
          
          if(i==1) c$i <- c$civil
          if(i==2) c$i <- ifelse(c$education >= 8,1,0 )
          if(i==3) c$i <- ifelse(c$gender == 1,1,0 )
          if(i==4) c$i <- ifelse(c$salary < c$med_salary,1,0 )
          if(i==5) c$i <- ifelse(c$year_hired >= 2005,1,0 )
          if(i==6) c$i <- ifelse(!is.na(c$old_party),1,0 )
          
          aa <- data.table(group_by(c,ibge,t,i) %>% summarise(outcome.new=sum(inter_rel),inter_tot=sum(inter_tot)))
          bb <- data.table(group_by(aa,ibge,i) %>% summarise(total.new=sum(inter_tot)))
          aa$inter_tot <- NULL
          
          newby <- left_join(v,aa,by=c("ibge","t","i"))
          newby[is.na(newby)] <- 0
          newby <- left_join(newby,bb,by=c("ibge","i"))
          newby[is.na(newby)] <- 0
          newby$var <- newby$outcome.new*100 / newby$total.new
          
          for (j in 0:1){
            wx <- filter(newby, i == j) 
           covars <- wx %>% select(age,gender,public_servant,newcomer,education,size,past_memb1,past_memb2,past_memb3,pt_coalition,pt_federal,pt_mayor)
            covars <- cbind(covars)
            
            reg <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars) 
            reg2 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars, level=90) 
       
            if(j==0){
            res[i,1] <- reg$coef[1]
            res[i,2:3] <- reg$ci[3,1:2] 
            res[i,4:5] <- reg2$ci[3,1:2] 
            }else{
            res2[i,1] <- reg$coef[1]
            res2[i,2:3] <- reg$ci[3,1:2] 
            res2[i,4:5] <- reg2$ci[3,1:2] 
              
            }
          }
        
        }
        
        
        # Do Plot
        res$type <- "a"
        res2$type <- "b"
        res$n <- res2$n <- 1:6
        res2$n <- res2$n+0.4
        
        resx <- rbind(res,res2)
        namesplot <- c("No Job Security","Job is Secure",
                       "Secondary Education","Post-Secondary Education",
                       "Male","Female",
                       "Low Salary","High Salary",
                       "Not Recently Hired","Recently Hired",
                       "Not Previously Partisan","Previously Partisan")
 
        resx <- resx %>% arrange(n)
        resx$nm <- namesplot
        
        ggplot(data=resx, aes(x=n, y=X1, group=type, colour=type))   + theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          geom_hline(yintercept=0, lty=3)+
           geom_point( aes(x=n, y=X1, group=type, colour=type), size=5)+
          scale_colour_manual(values=c("gray40","gray90"))+
          geom_errorbar(aes(ymin=X2,ymax=X3 , group=type, color=type),linewidth=1, alpha=.85, width=0) +
          geom_errorbar(aes(ymin=X4,ymax=X5 , group=type, color=type),linewidth=3, alpha=.65, width=0)+
          scale_x_continuous(breaks=unique(resx$n), labels=resx$nm)+
          theme(axis.text = element_text(colour="black", size=11, family="A")) +
          theme(axis.title = element_blank())+
          theme(legend.position = "none")+
          coord_flip()
        
         
        
        
         
        rm(c,v,wx,w1,w2,reg,res,resx,res2,namesplot,bands,b1,b2,covars,newby,aa,bb,reg2)
        gc()
        
        
    # =============================================================================================================================
         
             
    
        
      
     
    
       
   
      
     
 
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
        
      
  
      
      
      
      
    